---
title: "results of list experiment and balance check"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output: html_document
---

```{r setup, include=FALSE}
rm(list=ls())
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)

library(tidyverse)
library(haven)
library(stargazer)
library(BalanceR)
library(ggpubr)

dat <- read_dta("data/iraq2022_ver1.dta")
```

# data summary
```{r}
stargazer(as.data.frame(dat), type = "text")
```

# simple results
```{r}
prop.table(table(dat$violent_control))
prop.table(table(dat$violent_treatment))
prop.table(table(dat$violent_experiment))

prop.table(table(dat$PMU_control))
prop.table(table(dat$PMU_treatment))
prop.table(table(dat$PMU_experiment))
```

# balance violence
```{r}
balance_violence <- BalanceR(data = dat, group = list_violent,
                         cov = c(sex, age, education, income))
print(balance_violence, digits = 3)
summary(balance_violence, digits = 5)
plot(balance_violence, point.size = 5, text.size = 18)
```

# balance PMU
```{r}
balance_PMU <- BalanceR(data = dat, group = list_PMU,
                         cov = c(sex, age, education, income))
print(balance_PMU, digits = 3)
summary(balance_PMU, digits = 5)
plot(balance_PMU, point.size = 5, text.size = 18)
```

# list results
```{r}
mean(dat$violent_treatment, na.rm = TRUE) - mean(dat$violent_control, na.rm = TRUE)
mean(dat$PMU_treatment, na.rm = TRUE) - mean(dat$PMU_control, na.rm = TRUE)

t.test(dat$violent_control, dat$violent_treatment)
t.test(dat$PMU_control, dat$PMU_treatment)
```

# list results
```{r}
list_pmu <- dat %>% 
  summarise(control_mean = mean(PMU_control, na.rm = TRUE), control_sd = sd(PMU_control, na.rm = TRUE),
            treatment_mean = mean(PMU_treatment, na.rm = TRUE), treatment_sd = sd(PMU_treatment, na.rm = TRUE))

sd(dat$PMU_control, na.rm = TRUE) / sqrt(length(dat$PMU_control))
sd(dat$PMU_treatment, na.rm = TRUE) / sqrt(length(dat$PMU_treatment))

p1 <- data.frame(
  Experiment = c("control", "treatment"),
  Mean = c(1.901985, 2.313669),
  SE = c(0.01705567, 0.02435765)
)
p2 <- ggplot(p1, aes(x = Experiment, y = Mean, group = Experiment)) +
  geom_point(stat = "identity", color = "black") +
  geom_errorbar(aes(ymax = Mean+SE, ymin = Mean-SE), width = 0.2, position = position_dodge(.9)) +
  labs(title = "List expriment of PMU") +
  theme_bw()


list_vio <- dat %>% 
  summarise(control_mean = mean(violent_control, na.rm = TRUE), control_sd = sd(violent_control, na.rm = TRUE),
            treatment_mean = mean(violent_treatment, na.rm = TRUE), treatment_sd = sd(violent_treatment, na.rm = TRUE))

sd(dat$violent_control, na.rm = TRUE) / sqrt(length(dat$violent_control))
sd(dat$violent_treatment, na.rm = TRUE) / sqrt(length(dat$violent_treatment))


p3 <- data.frame(
  Experiment = c("control", "treatment"),
  Mean = c(1.730179, 2.199164),
  SE = c(0.01782132, 0.02394243)
)
p4 <- ggplot(p3, aes(x = Experiment, y = Mean, group = Experiment)) +
  geom_point(stat = "identity", color = "black") +
  geom_errorbar(aes(ymax = Mean+SE, ymin = Mean-SE), width = 0.2, position = position_dodge(.9)) +
  labs(title = "List expriment of election violent") +
  theme_bw()

p4
p2

ggarrange(p4, p2,
          nrow = 1, ncol = 2, align = "hv")
```
